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Abstract 

We compare the performance of the PHMC algorithm with the one of 
the HMC algorithm in practical simulations of lattice QCD. We show that 
the PHMC algorithm can lead to an acceleration of numerical simulations. 
It is demonstrated that the PHMC algorithm generates configurations car¬ 
rying small isolated eigenvalues of the lattice Dirac operator and hence 
leads to a sampling of configuration space that is different from that of the 
HMC algorithm. 
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Introduction 


In this paper we continue our discussion of the Polynomial Hybrid Monte Carlo 
(PHMC) algorithm [Q, This algorithm, designed for simulations of models 
containing fermionic degrees of freedom, is based on the idea 0 of combining 
the Hybrid Monte Carlo (HMC) algorithm ^ with the multiboson technique 
[^. In the PHMC algorithm the update part relies on an approximation of the 
exact fermion action to be simulated. The error induced by this approximation 
is corrected for by a reweighing technique, which introduces a correction factor 
taken into account in the sample average of the observables. 

In this paper we will present our results concerning the dynamical behaviour of 
the PHMC algorithm in practice. On the quantitative level we will compare its 
performance with the one of the HMC algorithm. Our numerical tests have been 
done in the Schrodinger functional set up p, |^, |[ , on small and moderately large 
physical volumes but at almost vanishing quark mass, which is feasible when 
using Schrodinger functional boundary conditions. We remark that since we are 
working at tiny values of the quark mass, the condition number of the fermion 
matrix employed in our simulations becomes 0(2000). 

We will demonstrate that the PHMC algorithm samples configuration space dif¬ 
ferently from the HMC algorithm. In particular, using the PHMC algorithm, 
gauge conhgurations with very small eigenvalues of the lattice Dirac operator can 
be reached. Consequences of this behaviour on the results for physical observables 
are discussed. 

We assume that the reader is familiar with refs. |l|, § . In particular, in the latter 
reference we have discussed a number of relevant technical aspects, which lay the 
ground for the present performance analysis. 

1 Numerical simulations with 
the PHMC algorithm 

In order to make the paper reasonably self-contained, we summarize here some 
features of the PHMC algorithm. We remark that throughout the paper we will 
use 0(a)-improved Wilson fermions. 
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1.1 Ingredients of the PHMC algorithm 

Denoting the lattice gauge link from x to x + a/t by U^{x) G SU (3) and a gauge 
field configuration by f/, the expectation value of any gauge invariant observable 
O = in full QCD with Uf = 2 degenerate flavours, may be written as 


(O) = z-^ 


I VUe-^^^^^det{Q^[U])0[U] , 


( 1 ) 


where Sg is the standard plaquette action for the pure gauge sector and Q is the 
Dirac operator for 0(a)-improved Wilson fermions multiplied by 75 (see below). 
The PHMC algorithm makes use of a polynomial approximation of The 

polynomial in a real variable s and having degree n is denoted by Pn,e{s) and 
constructed such that it approximates in the range 0 < e < s < 1 with a 
relative £t error bounded by 
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( 2 ) 


We choose the normalization of the Dirac operator such that the highest eigen¬ 
value of is smaller than 1 and write the corresponding polynomial in 
Pn,e{Q‘^)i in a factorized form: 


n 


k=l 


(3) 


where Cn,e is a positive constant, is determined by yTfc (see 0 ] for the exact 
relation), and the ZkS are the complex roots of Pn,e{,s). We note that special care 
has to be taken concerning the precise ordering of the factors in eq. (|]) in order 
to avoid problems with rounding errors |^. 

The full QCD (n/ = 2) partition function may now be represented as 

Z = j VUV(^^V(j)Vr]^Vr]WeP^^+^^+^^^^ 

Sp = 5p[[/,0] = (^tp„,,(Q2[[/])0 (4) 


by introducing the auxiliary pseudofermion fields (i.e. boson helds with spin and 
colour indices) 0, rj and the correction factor W = W[r], U]\ 

W = exp { 7^(1 - [Q^ ■ Pn,eiQ^)]~^)r]} ■ (5) 
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Each evaluation of W requires a trivial Gaussian “update” of the r^-held and 
the solution of the system [Q‘^Pn,^{Q‘^)]x = V- practice it turned out to be 
useful to generate the T^-helds Ncorr times for each given gauge held conhguration. 
Denoting averages evaluated with the effective action Sg + Sp + T]^r] as (...)p, the 
exact averages denoted as (...) are obtained by reweighing with W 

{O) = {W)-p\OW)p . (6) 


In |l[] we presented some tests of the PHMC algorithm for non-improved Wilson 
fermions. In this paper, we extend these tests to the case of 0(a)-improved 
actions. With respect to the non-improved case this amounts to adding the so- 
called “clover” term to the lattice Dirac operator, as specihed below. The 
modihcations in the PHMC algorithm induced by this extra term in the action 
are completely analogous to the ones needed for the standard HMC algorithm and 
our implementation of the PHMC algorithm for 0(a)-improved fermions follows 


closely the procedure described in for the HMC algorithm. 

For the actual simulations we consider hypercubic space-time lattices with lattice 
spacing a and size x T. With the lattice spacing set to unity from now on, 
the points on the lattice have integer coordinates (xq, xi,X 2 , X 3 ), which are in the 
range 0 < Xq < T-,0 < Xi < L. The gauge and the fermion helds obey Schrodinger 
functional boundary conditions as used in and detailed in [^, |^, |^. The matrix 
dehning the fermion action will be denoted by Q: 


Q{U),y = ^75|(1 + ElP. 


Cm 


J gyp 




^^{(1 + (1 + ; ( 7 ) 

where k is the hopping parameter and Cgw the improvement coefficient. The 
constant cm serves to optimize the simulation algorithm and cq = (1 For 

further unexplained notations we refer to refs. [|, 

In order to speed up the Monte Carlo simulation, not the original matrix Q but 
an even-odd preconditioned matrix Q is used. We expect the algorithm to 
be working equally well by using different preconditioning techniques like SSOR 
0 - Let us rewrite the matrix Q in eq. (0) as 



Meo 

1 + To( 


( 8 ) 
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( 9 ) 


where we introduce the matrix T^e {Too) on the even (odd) sites as 




The off-diagonal parts Meo and Moe connect the even with odd and the odd with 
even lattice sites, respectively. Preconditioning is now realized by writing the 
determinant of Q, apart from an irrelevant constant factor, as 


det(Q) oc det(l-|-Tee) det Q 


Co 


Q = ^75(1+ Toe -M,e(l+Tee)-Wee) • 


Cm 


( 10 ) 


The constant factor cq is given by Cq = 1/(1 + 64k^), and the constant cm is 
chosen such that the eigenvalues of Q are well within the interval [—1,1]. Since 
for the simulation algorithms the eigenvalues have to be positive, we dually work 
with the matrix Q^. We note that in the case Csw 7^ 0 the PHMC algorithm makes 
use of the polynomial approximation Pn^^iQ^) — 1/(0^) only to simulate detQ^, 
while the term det(l -|- Tee)^ is treated exactly. The correction factor therefore 
accounts only for the missing contribution, i.e. det Pn,e(^Q‘^), to the partition 
function. 


1.2 The simulations 

An important question is how the parameters of the polynomial Pn,e are to be 
chosen. Following ref. |^, a practical recipe for the choice of e and n may be 
given by 

^ ^ 2 (-^min(Q^)) ^2) 

(A„.ax(Q2)) 

and the value of n is set such that 5 ~ 0.01 (see eq. (|)). In eq. ([TT|) Amin(Q^) and 
Amax(Q^) denote the lowest and the highest eigenvalues of respectively. Our 
experience suggests that only a poor knowledge of the value of the average condi¬ 
tion number k = (Amax(Q^)/Amin(Q^)) for the specihc run parameters is needed. 
We remark that k ~ (Amax(Q^))/(Amin(Q^))- An estimate of k can be obtained, 
e.g. in the thermalization phase of the simulation, which may be performed by 
using either the standard HMC algorithm or the PHMC algorithm itself. We 
have also found that a very good and decisive check about the quality of the cho¬ 
sen polynomial approximation can be performed by monitoring the fluctuations 


5 



of the correction factor W: using too poor a polynomial approximation to Q~^ 
gives rise to large fluctuations of W, and consequently large fluctuations of many 
reweighted observables (eq. @), which can be detected after a few trajectories. 
Another remark [| concerns the dependence of the approximation on the volume 
V: the difference of the actions As = Shmc — *S'phmc is asymptotically As = 
V Cs exp(—2y^n), with Cs some proportionality constant. Since e is hxed by 
the condition of eq. (ED we find, if we also want to keep As fixed, that n ~ 
— (log A 5 — log Id — logC 5 )/ 2 -y/e. We see that the explicit volume dependence in 
n is rather weak in comparison with the (power-like) volume dependence induced 
by the way we choose e, following the criterion of eq. We therefore expect 

that the PHMC algorithm will also work efficiently in the case of large volumes, 
while keeping the value of 5, eq. (^, about 0 . 01 . 

The numerical tests are performed on 8^ x 16 lattices using the massively parallel 
Alenia Quadrics (APE) machines. Simulation parameters were chosen to be 

^ = 6.8 , K = 0.1343 , Csw = 1.4251 (12) 

^ = 5.4 , K = 0.1379 , Csw = 1.7275 . (13) 

These parameter values correspond to those used in simulations to determine the 
values of Cgw non-perturbatively 0 . 

All tests described below were performed on the APE machines by running Wep 
replica in parallel, with Wep set to 32 or 16. Since the Wep replica were indepen¬ 
dently thermalized, the data from the different replica are statistically indepen¬ 
dent from each other. This allows for a reliable error analysis, provided that for 
each replicum the statistics is several times larger than the integrated autocor¬ 
relation time of the observable considered. We determined our statistical errors 
for the observables, given below, from the variance of the Wep data obtained 
from running in parallel. We checked that the results were consistent with those 
obtained from a jack-knife procedure combined with binning. We refer to for 
further details. 

It is also possible to divide the Wep system replica into 2 sets of Wep/2 replica and 
analyse each of these two sets of data (a and b) separately. This gives two errors 
Aa and A;,, each of them obtained with half the statistics of the full run. By 
^We are grateful to A.D. Kennedy for this argument. 
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rescaling Aq and A;, by \/2, we can obtain an estimate of the error on the error, 
which in turn gives a measure of the error on the integrated autocorrelation time. 
This way of determining the error on the error yields values that are compatible 
with the estimate B0 of the relative error on the error given by (2Arep) 
or, in the case of a binning analysis leading to Abiock independent measurements, 
by (2iVbiock)~^^^- In the latter case, of course, the values of iVbiock have to be large 
enough that a plateau behaviour of the error can be detected. For a few tables 
below, besides the mean values and the errors (indicated in round brackets), we 
also quote the error on the error (indicated in square brackets). 

2 Results at /? = 6.8 

In this section we give our results for various quantities as computed with the 
PHMC algorithm and compare with those obtained from the HMC algorithm. We 
will compare bulk quantities as well as quark correlation functions and certain 
combinations of them. 

We give in table |I| the parameters for both simulation algorithms. As reported 
in [0, in the simulations with the HMC algorithm, we found that sometimes a 
trajectory was not accepted for a number of times. The cure was that every l-th 
trajectory the step size 6t was changed to a smaller value, and the corresponding 
number of molecular dynamics steps, Amd, was increased to reach a unit trajec¬ 
tory length. In the actual simulation a value of / = 6 was chosen and we give in 
table 1^ the effective values of 6 t and Amd built from the normal and the smaller 
step size. We remark that this effect, observed in simulations with the HMC al¬ 
gorithm, never appeared within the simulations using the PHMC algorithm and 
that the step size was always kept constant there. This allowed in particular to 
run the PHMC algorithm at an acceptance rate smaller than the one obtained 
with the HMC algorithm. 

2.1 Bulk quantities 

As bulk quantities we consider the expectation values for the plaquette P, the 
lowest Amin and the largest Amax eigenvalues of Q^, and the derivative of the pure 
gauge action with respect to the background held, dSg/drj. The latter quantity 
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Table 1: The parameters for both simulation algorithms aX (3 = 6.8. We denote 
by Stat the statistics obtained in units of trajectories and Pace is the acceptance 
rate. 


Algorithm 



P 

^ acc 

Stat 

e 

n 

Cm 

HMC 

0.059 

16.6 

0.948(8) 

2944 

- 

- 

- 

PHMC 

0.077 

13 

0.79(1) 

2944 

0.0022 

62 

0.735 

PHMC* 

0.077 

13 

0.758(8) 

2944 

0.0022 

54 

0.725 


Table 2: Comparison of bulk quantities as obtained from the two algorithms. We 
use the notation PHMC(iVcorr) to indicate the values of Worr used for the analysis. 
Ncorr = 0 means that the correction factor is set to 1. In square brackets we give 
our estimate of the error on the error. 


Algorithm (P) 

HMC 0.673384(53) [7] 

PHMC(4) 0.673483(46) [6] 

PHMC(2) 0.673496(45) [6] 

PHMC(l) 0.673505(48) [6] 

PHMC(O) 0.673512(44) [6] 

PHMC*(2) 0.673435(66) [8] 


(An,in(Q')) (A^ax(Q')) {dSJdr,) 


0.001150(35) [4] 
0.001152(42) [5] 
0.001150(43) [5] 
0.001141(42) [5] 
0.001025(46) [6] 


0.87188(25) [3] 
0.87164(36) [5] 
0.87173(39) [5] 
0.87190(42) [5] 
0.87177(30) [4] 


23.1(2.4)[0.3] 

22.1(2.3)[0.3] 

22.6(2.3)[0.3] 

22.6(2.2)[0.3] 

20.6(2.0)[0.3] 


0.001117(44) [6] 0.87426(70) [9] 27.1(3.1) [0.4] 


can be used to dehne a running coupling constant in the pure gauge theory 
As table shows, we hnd that, for Worr > 0 the values of all bulk quantities 
are completely consistent with the corresponding ones from the HMC run. Also 
the uncorrected (see Worr = 0) values for (Amax) and {dSg/drj) are in agreement 
with the HMC values while, perhaps, the ones for (P) and Amin are somewhat 
off. Within the error on the error, also the estimated errors on the observables 
are consistent between the PHMC and HMC algorithm. 

The prominent exception is Amax, where the error from the PHMC algorithm 
appears to be substantially larger than the one from the HMC algorithm. Note, 
however, that the mean value and the error for the uncorrected value of Amax are 
both consistent with the corresponding quantities from the HMC algorithm. In 
addition, the error decreases when Worr is increased from 1 to 4. This points 
















towards the interpretation that the larger error is jnst indnced by the additional 
noise appearing throngh the correction factor and that there is no large antocor- 
relation time hidden in the PHMC algorithm. Of course, Amax is a pure cut-off 
quantity and is not expected to be physically relevant. 


2.2 Quark correlation functions 

Quark correlation functions are important quantities, from which many physical 
observables can be constructed. We hence extend our comparison of the algo¬ 
rithms by considering certain quark correlation functions, which are often used 
in computations with Schrodinger functional boundary conditions. To this end 
we closely follow ref. and construct correlation functions using boundary quark 
helds (, ( at Euclidean time xq = 0: 

fA{xo) = -^iAl“(x)C(0,y)75^r“C(0,z) 

y,z ^ ^ 

fp{xo) = -^^P“(a;)C(0,y)75^r“C(0,z) . (14) 

y,z ^ 

In eq. (p!4D Aq{x) denotes the isovector axial current and P°‘{x) the corresponding 
density 

= V'75^r>, (15) 

where r“ is a Pauli matrix acting on the flavour indices of the quark held. 
Analogously one may build f'j^{T — xq) and f'p{T — xq) with boundary quark helds 
C', C cl-t Xq = T. 

We will consider the correlation functions fp{xo) as well as hnite diher- 

ences of them: 


dA{xo) = {Oq -F do)fA{xo) , 0 < xo < T 
Dp{xo) = dQdofp{xo) , 0 < xo < T . 


(16) 


In eq. (|^) do is the lattice forward derivative, and the lattice backward deriva¬ 
tive 


= f{xo + 1 ) - f{xo) 

= f{xo) - /(xo - 1) . 


dof{xo) 

dofixo) 


9 


( 17 ) 



Table 3: Comparison of quark correlation functions as obtained from the two 


algorithms. Notations are as in table 

Algorithm 

(/^(T/2)) ■ 10-4 

(/p(T/2)) ■ 10-6 

(d^(T/2)) ■ 10-4 

^ (Dp(T/2)) ■ 10-6 

HMC 

0.542(39) [5] 

0.1072(55) [7] 

-0.171(13)[2] 

0.2062(60) [8] 

PHMC(4) 

0.609(58) [7] 

0.1074(51) [6] 

-0.171(12)[2] 

0.1957(55) [7] 

PHMC(2) 

0.612(60) [8] 

0.1075(53) [7] 

-0.168(13)[2] 

0.1969(56) [7] 

PHMC(l) 

0.618(60) [8] 

0.1081(53) [7] 

-0.170(13)[2] 

0.1986(56) [7] 

PHMC(O) 

0.879(144) [18] 

0.1300(80) [10] 

-0.193(15)[2] 

0.1891(83)[10] 

PHMC* (2) 

0.632(70) [9] 

0.1088(56) [7] 

-0.173(12)[2] 

0.1939(50) [6] 


We compare the results for quark correlation functions as obtained from the two 
algorithms in table |^. We take the distance in time to be half the temporal size of 
the lattice, i.e. Xq = T/2. We can see from table that we have again consistent 
results for the mean values of /p(T/2), d^(T/2) and Dp{T/2), as well as for the 
corresponding errors. However, for /a(T/ 2) the error from the PHMC algorithm 
is substantially larger than the one from the HMC algorithm. The discrepancy 
is well outside the uncertainty of the error as indicated by the error on the error 
given in the square brackets. Even more pronounced is the behaviour of the 
uncorrected value of /a(T/ 2). The error is a factor of about 4 larger than in the 
HMC case, and also the mean value is off. 

A hrst step for understanding the larger error obtained from the PHMC algorithm 
is to look at the distribution of /yi(T/2), which we show in Eg. [^. It is clearly 
seen that the distribution as obtained with the PHMC algorithm spreads much 
further out, towards large values of /^(T/2). In principle, this effect can come 
either from a large autocorrelation time encountered within the PHMC algorithm 
or from a different sampling in conhguration space. To decide on this question, 
we plot in Eg. ^ the time evolutions of various quantities. 

Let us start with the correction factor W. Figure ^(a) shows that W can be¬ 
come small, assuming values that are clearly much below the average value, 
(W) ~ 0.45. At the time when IT -C 1, /a(T/ 2) assumes very large values 
as shown in Eg. ^b). In Eg. ^(c) we show W ■ /a(T/ 2): the spike in Ja is now 
suppressed by the correction factor W. The PHMC algorithm seems to allow for 
conEgurations that make large contributions to quark correlation functions and 
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0 2x10“ 4x10“ 6x10“ 

f,(T/2) 


Figure 1: The distributions of the quark correlation function /^(T/2) 
as obtained from the HMC and the PHMC algorithms at /d = 6.8. 
Note that for the PHMC algorithm we plot the uncorrected values for 
/n(T/2). 


partly suppresses these contributions in the reweighing procedure through small 
values of the (noisy) correction factor. In £g. §(d) we show the Monte Carlo time 
evolution of /^(T/2) for the HMC algorithm, which looks quite different from 
that of the PHMC algorithm. We conclude that the difference in the variance of 
/yi(r/2) is not due to a large autocorrelation time but reflects the fact that the 
PHMC algorithm really samples the conhguration space differently. A similar 
observation was made in in a different context. 

The Monte Carlo time evolution of /p(T/2) is plotted in hgs. ||(c,d) for the 
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Figure 2: The Monte Carlo time evolutions for W, /^(T/2), 

WfA{T/2), as obtained from the PHMC algorithm, and /^(T/2) from 
the HMC algorithm. 
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Figure 3: The Monte Carlo time evolutions for the correction factor 
hF, Amin(O^)) /p(^/2), as obtained from the PHMC algorithm, and 
fp{T/2) from the HMC algorithm. 
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PHMC and the HMC algorithms. For this quantity we do not observe spikes in 
the PHMC Monte Carlo history and both time evolutions are similar. This is 
consistent with the results in table 1, where we saw that the errors for /p(T/2) 
are comparable in the two cases. In figs. |5(b,a) we show the time evolution of the 
lowest eigenvalue Aniin(Q^) and W for the PHMC algorithm. 


2.3 Combinations of quark correlation functions 


Following H, 121 we dehne correlation functions 


r(xo) = ^(do +do)fA(xo)/fp(xo) 
= ^dodofp(xo)/fp(xo) 


(18) 


and analogously r'(T — xq), s'(T — Xq) in terms of f'^iT — Xq) and fp(T — xq). 
These correlation functions allow us to dehne an unrenormalized PCAC current 
quark mass: 

M{xo, yo) = r(xo) - (19) 

s{yo) - s'{yo) 

and analogously M'. The non-vanishing of the difference between M and M' at 
certain time distances 


AM = M(jr, ir) - M'( jT. jT) (20) 

is a lattice artefact appearing linear in the lattice spacing. The requirement that 
AM assumes its tree-level value, AM = 0.000277, is the improvement condition 
to determine the values of Cgw non-perturbatively. 

We may build various, physically interesting combinations of the correlation func¬ 
tions of eqs. (^). We will consider the unrenormalized current quark mass M 
(eq. (1^), AM (eq. (^0|) ), and an estimator oi the improvement coefficient ca, 

. r(r/4) - r'(T/i) 

s(T/4) - s'(T/4) ■ ' ’ 

We want to emphasize that ca should not be taken as the true non-perturbatively 
determined values of ca- We consider ca in this work as a purely technical 
parameter, which can also be used in comparing the two algorithms. We give 
our results for M, AM and ca in table We hnd that, at least within the 
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Table 4: Comparison of combinations of quark correlation functions as obtained 
from both algorithms. Notations are as in table ^ 


Algorithm ca M 


HMC 

PHMC(4) 

PHMC(2) 

PHMC(l) 

PHMC(O) 

PHMC*(2) 


-0.0265(28) [4] 
-0.0262(30) [4] 
-0.0257(26) [3] 
-0.0265(31) [4] 
-0.0242(27) [3] 
-0.0247(56) [7] 


0.00144(36) [5] 
0.00161(35) [4] 
0.00155(33) [4] 
0.00150(40) [5] 
0.00194(27) [3] 
0.00180(62) [8] 


AM ■ 103 

0.045(311) [40] 
-0.086(390) [50] 
-0.129(366) [50] 
0.015(381) [50] 
-0.020(493) [60] 
0.745(404) [50] 


statistical uncertainties, the average values as well as the errors are completely 
compatible. 

We close this section by remarking that we also performed a simulation with the 
PHMC algorithm choosing a trajectory length of Amd^'r ~ 0.5. The results from 
this test are, however, rather inconclusive: while for some observables the errors 
did not change with respect to the run with unit trajectory, for other observables 
we found an increase of the errors as expected. 

3 Results at /3 = 5.4 

This section is devoted to a discussion of the results obtained at /? = 5.4, for 
which the lattice spacing a ~ 0.1 fm. We set k = 0.1379 and Cs* = 1.7275. At 
these values of the parameters we hnd a quark mass M = 0.009(1) We use 
a 8^ X 16 lattice and the boundary conditions are the same as in section 2. For 
reasons that will become clear from our discussion, we do not aim in this section 
at a comparison of the HMC and PHMC algorithms on the same quantitative 
level as it was done in the previous section for the results at /3 = 6.8. We will 
rather emphasize the qualitative behaviour of the PHMC algorithm in sampling 
conhguration space and reweighing observables when very small eigenvalues of 
occur. 
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3.1 Low-lying eigenvalues 


As was shown in p| for the parameter values considered in this section, isolated 
very small eigenvalues of the operator are found. We illustrate this again in 
£g. 1^, by showing the Monte Carlo time evolution of the five lowest eigenvalues in 
four typical situations. In £g. ^(a) the hve lowest eigenvalues lie in a narrow band 
and we hnd a basically continuous spectrum, at least up to the tenth eigenvalue. 
In figs.^(b,c) there are a few eigenvalues that assume rather small values and 
hnally, in £g. §(d), we observe very small, isolated eigenvalues, lying many orders 
of magnitude below the ones in £g. H(a). As can be seen from the distribution of 
Amin in £g. 3 of ref. p|, such very small eigenvalues could not be observed in the 
corresponding simulations using the HMC algorithm. 

It is a natural question to ask, whether the occurrence of the small eigenvalues 
shown in £g. ^ is related to some topological effects. We therefore consider the 
values of the pure gauge action and the naive topological charge after per¬ 
forming 500 cooling iterations (see also [|^]); these values will be denoted in 
the following by S'dassicai and Qtopo, respectively. We emphasize that we do not 
want to give a precise and reliable number for the topological charge itself, but 
rather that we are interested in the qualitative behaviour of Qtopo and in only 
estimating the autocorrelation time of a quantity that is related to topology. We 
remark that in the case of Schrodinger functional boundary conditions there exist 
bounds P on the pure gauge action Sg given by 


goSg > vr^ Qtopo = 0 , 

g^Sg > Svr^lQtopol • 


( 22 ) 


In £g. D we plot an example of the Monte Carlo time evolution of S'dassicah Qtopo 
and the lowest eigenvalue of Q^. It is remarkable that, although working at 
basically zero quark mass, we see some transitions between topological sectors. 
As expected from the bounds of eq. (^^ the behaviour of S'dassicai closely follows 


the one of Qtopo- 

The behaviour of the lowest eigenvalue of and the topological charge are not as 
closely related. Small eigenvalues are expected when a transition between topo¬ 
logical sectors occur and, indeed, there is one instance, shown in £g. |^, where 
exactly this happens. However, we also see from £g. that the topological charge 
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Figure 4: We show the Monte Carlo time evolution of the hve lowest 
eigenvalues of at /9 = 5.4 in four typical situations. The lowest 
eigenvalue is shown by the open symbols, the remaining eigenvalues 
by the hlled ones. 
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can change withont any occurrence of a very small eigenvalue. This might of 
course be due to the fact that a small eigenvalue has appeared during the mo¬ 
lecular dynamics evolution. Finally we can have situations where the eigenvalue 
becomes very small but the topological charge does not change, which might 
correspond just to an unsuccessful attempt to change topological sectors. The 
relation between the topological charge and very small eigenvalues may be partly 
obscured in our case by the fact that we use only a naive definition of the topo¬ 
logical charge. Moreover we remark that the index theorem does not have to 
hold owing to the existence of lattice artefacts and to our choice of Schrodinger 
functional boundary conditions. In any case, since we are working at almost zero 
quark mass and reasonably large physical volume, we take fig. ^ as an encourag¬ 
ing indication that the PHMC algorithm is able -even in this physical situation- 
to explore different topological sectors. Of course, only more extensive investi¬ 
gations, possibly at larger physical volumes, can decide whether our tentative 
conclusion is too optimistic. 

We remark that when measuring the topological charge with the PHMC algo¬ 
rithm, its physical value will be the one reweighted with the correction factor. If 
we are close enough to the continuum for the effects of the lattice spacing to be 
negligible and we are able to work at vanishing quark mass, a non-trivial topo¬ 
logical charge has to induce the appearance of a zero mode. Since the correction 
factor is proportional to this zero eigenvalue, the reweighted topological charge 
will always be zero. This corresponds, of course, exactly to the continuum situ¬ 
ation, where the topological charge vanishes after integration over the fermions, 
provided that at least one of the fermion species is massless. 

3.2 Modified correction factor 

As discussed in [0, in order to deal with situations where very small eigenvalues 
may occur, the correction factor W, eq. (|), has to be modified. The reason is that 
in the presence of very small eigenvalues the noisy estimate of det[Q‘^Pn,e{Q^)] 
given in eq. (^) is largely dominated by those p-fields that are almost orthogonal to 
all the eigenfunctions of the small eigenvalues. Since the probability of extracting 
such //-fields from a distribution cx exp (—//I//) is low, we would need a large value 
of Ncorr to obtain a good (i.e. not too noisy) estimate of det[Q^Pn^e{Q^)]- 
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An improved definition of the correction factor, replacing eq. (|^), is given by (see 
P for more details): 

W = fVBhhm . (23) 

The separation between Wb and Win is controlled by a new parameter: e <C e. In 
eq. ( [23| ) Wb is a “bulk” factor, taking into account the contribution of all modes 
with eigenvalues larger than e: 

Wb[v, U] = exp {vi[Rn,e ■ {Q^ ■ Pn,e)-^]v±} (24) 

and WiB an “infrared” factor that incorporates the contribution from the eigen- 
modes of lying below e, 

W^IR= n [1 + ^n,.(A,)] . (25) 

Xj<e 

In eqs. (P^) and (P5|) we have introduced the relative fit deviation = Q^Pn,e — 
1, the eigenmodes |Aj) of and the projection of the //-held onto the subspace 
orthogonal to all the modes lying below e: 

Q2|A,) = A,|A,) 

\V±) = 1^)Aj)|Aj)(Aj|77) . (26) 

j 

Whereas Wb is given again by a noisy estimator, ITir is calculated “exactly” 
in terms of the eigenvalues of that are smaller than e. These eigenvalues 
can be explicitly computed, together with the corresponding eigenvectors, with a 
pre-dehned accuracy |]T^ 1^ order to guarantee the exactness of the PHMC 

algorithm, e has to be fixed in a simulation beforehand. For the present inves¬ 
tigation we have chosen e = e/10. Clearly, when no eigenvalues smaller than e 
occur, WiB = 1. In particular, for e = 0 we are back to the old correction factor 
of eq. (1). 

The difference between the old and the modihed correction factors, as evaluated 
on a gauge conhguration carrying a very small isolated eigenvalue (Amin = 3.7 ■ 
10“^), is demonstrated in fig. p. There we plot the distribution for a fixed gauge 
held conhguration of tc = log(IFBlFiR) as obtained from the generation of 600 
? 7 -helds for two diherent values of e. When setting e = 0 (open squares) the 
distribution is very broad, leading to a very noisy and imprecise estimate of 
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det[Q^-Pn,e(Q^)]- When setting e = 0.00011 (filled sqnares), i.e. a value ten times 
smaller than e, the distribution appears as a needle on the scale of the upper plot 
in fig. p. In the lower plot of this figure we therefore resolve the distribution for 
e = 0.00011. The picture nicely demonstrates that if we use the original form 
of the correction factor, i.e. set e = 0, the estimate of W[r]^U] is dominated 
by the terms W-Q[rjm,U] with (r^mlAmin) — 0. However, vectors \ri), which are 
almost orthogonal to |Amin), may be extracted very rarely from the probability 
distribution P[7]] = exp{—This leads to large fluctuations affecting the 
estimate of det[Q^Pn^e{Q‘^)], and a very large value of Worr is needed for the 
result to be sufficiently precise. For e = 0.00011 the distribution becomes much 
narrower and a value of Worr lower than 10 is sufficient to achieve a precision 
that is appropriate for the purpose of keeping the fiuctuations of Wb small. 

For situations where no eigenvalue of is exceptionally small it should make no 
difference whether e is set to zero or to some finite value smaller than e. The noise 
in the estimate of det [Q^Pn^eiQ^)] will be essentially the same for both cases, since 
there is no single mode that plays a dominating role in determining the value of 
WirWb- We have checked this expectation explicitly and our numerical results 
fully confirm the above picture. We have also checked that a relative precision of 
1% in the evaluation of the low-lying eigenvalues of yields eigenvectors that are 
accurate enough to get a precision sufficient for the projection onto the subspace 
orthogonal to the one spanned by the eigenvectors themselves. Concerning the 


computational cost of the modified correction factor, eq. (|^) , an overhead with 
respect to the cost of computing the ordinary correction factor, eq. (^), comes 


from the evaluation of the needed eigenvalues and eigenvectors of Q^. This over¬ 
head depends on the choice of e. In our test run at /3 = 5.4, we found that for 
a case (see below) when the four lowest eigenvalues and eigenvectors of are 
needed, the overhead for the modified correction factor is just half the time of 
evaluating the ordinary correction factor having Worr = 4. We mention that, 
when setting e = e/10, the modified correction factor had only to be computed 
in about 35% of our measurements. This leads to an additional reduction of the 
overhead. We will hence neglect this overhead when discussing computational 
costs in section 4. 

In table § we show data for the low end of the spectrum of Q^\ for the ten 
lowest eigenvalues, we consider the expectation values and the variance of the 
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Figure 6; The distribution of w = logiWBWin) for e = 0 (open 
squares) and e = 0.00011 (filled squares), on a fixed gauge configura¬ 
tion carrying an exceptionally small (isolated) eigenvalue of Q^. In the 
lower hgure, we resolve the distribution of w for the case e = 0.00011. 
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Table 5: The uncorrected (iVcorr = 0) values of the ten lowest eigenvalues of Q^: 
We give the expectation values, with the corresponding true error in parenthesis, 
and the variance, as obtained from our PHMC test at /3 = 5.4, see table Note 
that the value of e was set to 0.0011. Moreover we show the ten lowest eigenvalues 
of for two particular gauge configurations (Ci and C 2 ), the first of which has 
a very small value of Amin- 


Eigenvalue 

(A) 

t/(A^) - (A)^ 

Cl 

C2 

Al = Amin 

0.00032(5) 

0.00024 

0.0000017 

0.00052 

A2 

0.00054(5) 

0.00026 

0.00027 

0.00087 

A3 

0.00090(5) 

0.00034 

0.00041 

0.00137 

A4 

0.00114(6) 

0.00032 

0.00077 

0.00152 

A5 

0.00140(6) 

0.00033 

0.00098 

0.00170 

Ae 

0.00162(6) 

0.00033 

0.00139 

0.00171 

Ay 

0.00190(6) 

0.00032 

0.00141 

0.00189 

-^8 

0.00212(5) 

0.00032 

0.00204 

0.00217 

Ag 

0.00237(5) 

0.00031 

0.00206 

0.00262 

Aio 

0.00256(5) 

0.00031 

0.00260 

0.00274 


uncorrected (Worr = 0) eigenvalues. We see that the variance is almost constant 
and takes a value of the same order of magnitude as the average lowest eigenvalue 
of Q^. We also give the example of two particular gauge configurations, one with 
an exceptionally small eigenvalue and another with no exceptional eigenvalues. 
Note that for the hrst conhguration (Ci) all the eigenvalues Xj, with 1 < j < 
10, lie somewhat below the corresponding eigenvalues measured for the second 
conhguration (C 2 ). We infer from the results for the variance that for practically 
all gauge conhgurations of our sample there are only very few eigenvalues lying 
below e. This also justihes our choice of e = e/10 for the modihed correction 
factor. We remark that with this choice of e for evaluating the modihed correction 
factor, eq. (1^), we need not more than the four lowest modes of (see table H). 
Let us hnally demonstrate that, despite the diherent behaviour of the two algo¬ 
rithms in sampling conhguration space, compatible results are found within the 
present statistical uncertainties. In table ^ we give the algorithmic parameters 
for the simulations performed at /? = 5.4 as well as the acceptance rates and the 
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statistics. 


Table 6: Technical parameters for the algorithms aX (3 = 5.4. 


Algorithm 

5t 

Nrad 

P 

acc 

Stat e n Cm 

HMC 

0.032 

34 

0.948(8) 

5120 - - - 

PHMC 

0.056 

18 

0.83(1) 

1632 0.0011 76 0.806 


In table we present a comparison of the bnlk qnantities. Note that the statistics 
for the HMC rnn is abont a factor of 3 larger. We emphasize again, however, 
that with this small statistics the error on the error is snbstantial and no real 
comparison of the performance between the two algorithms is possible. To really 
say something abont the performance, a mnch larger statistical sample wonld be 
necessary for both algorithms. Since a non-negligible amonnt of compnter time 
has already been invested in obtaining the present statistics, we feel that snch 
a comparison shonld be made within a project that aims at the same time at 
physical resnlts. 

Table 7: Comparison of bnlk qnantities as obtained from the HMC and the PHMC 
algorithms at /3 = 5.4. Notations are as in table ^ 


Algorithm 

(P) 

(A„.in(Q')) 

(An,ax(Q')) 

{dSJdr,) 

HMC 

0.563331(65) 

0.000561(17) 

0.83555(31) 

0.8(2.0) 

PHMC(8) 

0.563302(120) 

0.000506(50) 

0.83672(99) 

-6.2(4.4) 

PHMC(4) 

0.563344(135) 

0.000528(69) 

0.83665(90) 

-3.6(5.0) 

PHMC(2) 

0.563404(138) 

0.000554(85) 

0.83649(190) 

-4.6(6.3) 

PHMC(l) 

0.563679(377) 

0.000600(107) 

0.83599(259) 

-0.7(9.6) 

PHMC(O) 

0.563336(122) 

0.000322(50) 

0.83730(43) 

-6.9(3.4) 


We remark that we have also monitored the qnark correlation fnnctions, intro- 
dnced in sections 2.2 and 2.3, hnding consistent resnlts for the HMC and the 
PHMC algorithm. No new qnalitative features arise with respect to our discus¬ 
sion for the data obtained at /3 = 6.8 (see the previous section). In particular, we 
hnd again spikes in the uncorrected quark correlation functions, in coincidence 
with gauge conhgurations carrying very small eigenvalues of Q^. For these con- 
hgurations we observe e.g. for /yi(T/2) values up to three orders of magnitudes 
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larger than the typical values assumed for “normal” gauge conhgurations; at the 
same time the modified correction factor, eq. (^), takes values up to three orders 
of magnitudes smaller than the usual ones, leading, as expected, to contributions 
of “normal” size to the reweighted average, eq. (|]). 

Taking A^corr = 8 and a statistics of 1632 trajectories for the PHMC algorithm, 
we find for the quark mass M = 0.0066(21) and for the lattice artefact AM = 
0.00299(183). With the same statistics, the HMC algorithm gives M = 0.0086(28) 
and AM = 0.00026(201). This indicates, but does not prove, that with the same 
statistics compatible errors can be obtained from the two algorithms also for these 
quantities. 

4 Computational cost 

A crucial question is, of course, how the cost of the PHMC algorithm compares 
with the one of the HMC algorithm. In this section we will therefore give the 
computational cost of both algorithms for generating one gauge held conhguration 
at the two values of (3 considered in this paper. For the simulations performed at 
(3 = 6.8, this comparison of the cost corresponds to a comparison of the actual 
cost to generate an independent conhguration, because the errors on almost all 
observables are compatible between the two algorithms when the same statistics 
is employed. For the simulation performed at /3 = 5.4, the situation is, however, 
diherent since, with the available statistics, the uncertainties on the integrated 
autocorrelation times are rather large and no dehnite statement can be made. 
However, regarding observables for which the very small eigenvalues are impor¬ 
tant, a comparison of the errors would be difficult even if the statistics were large. 
If the modes corresponding to these very small eigenvalues are physically impor¬ 
tant for some observables and the HMC algorithm generates these modes only 
very seldom, a direct comparison of the huctuations of these particular observ¬ 
ables computed with the two algorithms is not appropriate. This is, of course, 
a general problem when comparing algorithms with different behaviour in sam¬ 
pling conhguration space 0. In such a situation the algorithms have very different 

^One example would be the behaviour of the cluster and the Metropolis algorithms at a 
first-order phase transition. 
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autocorrelation times. 

In 0 we gave a detailed description of the computational cost of the PHMC 
algorithm in units of matrix times vector Q(j) operations. Therefore, we list here 
only the formulae for the cost analysis derived in . Let us remark that the cost 
of a single trajectory in both algorithms may be written as 

Ctot = Cq^ + Cextra , (27) 

where the hrst contribution is given by the number of matrix times vector Qcj) 
operations and the second part accounts for all other operations. Asymptotically, 
when the condition number of Q becomes large, Cq^ will by far dominate the cost 
of the algorithms. We will therefore only discuss and compare the cost in 
the following. 

Let us denote by Ncg fhe average number of iterations of the Conjugate Gradient 
algorithm that is implemented in our programs for all matrix inversions^. Then 
the cost for the HMC algorithm in units of Q(j) operations is given by 

Gq^(HMC) = 2 ■ {2N^^ + 1) ■ Ncg , (28) 

The factor (2Amd + 1) originates from the use of the Sexton-Weingarten inte¬ 
gration scheme |^. The cost for the PHMC algorithm is split into three parts 

E- 

Co*(PHMC) — Gbhb T Gupdate T Gcorr ) (29) 

where Gbhb is the cost for the heatbath of the bosonic held 0, Gupdate the cost for 
the computation of the variation of the action with respect to the gauge held and 
Gcorr the cost to evaluate the correction factor. In units of Q0 operations we hnd 

Gbhb = (2n + 2)-iV^g’ + n 
Gupdate = 3n ■ (2Amd + 1) 

Gcorr = (2n + 2) ■ ■ iVcorr . (30) 

The factor Ncorr denotes as usual the number of evaluations of the correction 
factor W per full gauge held update (or molecular dynamics trajectory). We 

^We remark that in the set up we consider here and using APE computers the standard CG 
solver was found to be competitive with the BiCGStab solver for the purpose of inverting Q^. 
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explicitly verified that the cost in real time, as expected from our formulae for 
Cupdate, C'bhb and Ccorr, agree with the one measured for our implementation of 
the PHMC algorithm on the APE computer. 


Table 8: Conjugate Gradient iterations and degree of the polynomial used in the 
PHMC runs. Notations are explained in the text. 



HMG 

PHMG 

/5 

NEg 

Ncg 

n 

iVbhb 

N^g 

6.8 

149.0(1) 

113.3(4) 

62 

3.66(4) 

3.26(4) 

6.8 

- 

- 

54 

3.61(6) 

3.88(3) 

5.4 

197.6(1) 

143.0(8) 

76 

3.56(6) 

3.88(4) 


All of the simulations done at /3 = 6.8 and [3 = 5.4 have been performed by 
running several replica in parallel. In particular for the HMC runs we always 
had 32 replica. Because the APE computer we are using is a SIMD machine, all 
replica have to wait until the Conjugate Gradient solver of the slowest replicum 
has converged. This “parallelization effect” has an important consequence for the 
HMC algorithm. We give in table the maximal number of CG iterations, 
as determined from the slowest replicum and the number of CG iterations Nqq, 
obtained by averaging over all replica. As we see from the tables, in particular 
for (3 = 5.4, there can be a substantial increase of the number of GG iterations 
from this parallelization effect. The analogous effect is much less relevant in the 
case of the PHMG algorithm, since it may occur only in C'bhb and Ccorr, which are 
asymptotically marginal in comparison with Cupdate- To be conservative in the 
estimate of the computational cost for the PHMG algorithm, we will neglect to 
correct for this small parallelization effect. We do mention, however, that doing 
so may reduce the values for Cbhb and Ccorr by a factor of 2 at /? = 5.4. 

From tables 0, I and 1^ we can now calculate the computational cost for both 
algorithms. We present the results in table ^ for (3 = 6.8 and in table |T^ for 
(3 = 5.4. We give the global costs for both algorithms considering the case 
of 32 replica (Cq^), where the HMG algorithm is slowed down by a significant 
parallelization effect, and the case of a single lattice system (Cg,^). 

For (3 = 6.8 we see that the dominating effect in the cost gain of the PHMG al¬ 
gorithm stems from the parallelization effect. Taking this effect out, we still have 
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Table 9: Computational cost for (3 = 6.8. We take Worr = 1 for the PHMC run 
with n = 62 (PHMC) and Worr = 2 for the PHMC run with n = 54 (PHMC*). 
The cost takes the parallelization effect into account when running 32 replica 
in parallel. Cq^ would be the cost when simulating a single lattice system. 


Algorithm 

Cbhb 

C^update 

r 

^corr 

^32 

Cq^ 

HMC 

— 

10192 

— 

10192 

7750 

PHMC 

523 

5022 

411 

5956 

5956 

PHMC* 

451 

4374 

854 

5679 

5679 

Table 10: Pure computational cost for [3 = 5.4. We consider the cases Worr = 4 
(PHMC(4)) and Worr = 8 (PHMC(8)). Notations are as in table ||. 

Algorithm 

Cbhb 

C^update 

r 

^corr 

/^32 

Cq^ 

HMC 

— 

27269 

— 

27269 

19734 

PHMC(4) 

624 

8436 

2390 

11450 

11450 

PHMC(8) 

624 

8436 

4780 

13840 

13840 


a performance of the PHMC algorithm better than that of the HMC algorithm, 
but the gain becomes marginal. We remark that at /3 = 6.8 the lattice spacing 
is very small and we are hence working in a correspondingly small physical vol¬ 
ume. Going to a more challenging situation, i.e. (3 = 5.4, we still hnd a large 
parallelization effect but now even if this is taken out, a factor of almost 2 is 
found in favour of the PHMC algorithm. We emphasize again at this point that 
we give here only the computational cost of the algorithms and do not take the 
autocorrelation time into account for the reasons discussed above. 


5 Conclusions 

In this paper we have tested the PHMC algorithm for 0(a)-improved Wilson 
fermions. We compared the computational cost of the PHMC algorithm, as 
well as its qualitative behaviour, with those of the HMC algorithm. Practical 
simulations were performed on 8^ x 16 lattices at /? = 6.8, which corresponds 
to a very small physical volume, and (3 = 5.4, corresponding to an intermediate 


28 










physical volume, with a lattice spacing a ~ 0.1 fm. The results of our tests lead 
us to the following conclusions: 

1) It is easy to hud values for the degree n and the infrared parameter e, deter¬ 
mining the polynomial approximation used in the PHMC algorithm, such 
that its performance becomes comparable to that of the HMC algorithm. 
As a guideline one may choose e ~ 2(Amin), with Amin the lowest eigen¬ 
value of the fermion matrix used in the simulation. The degree n of the 
polynomial should then be chosen such that 5 < 0.01, see eq. (^. 

2) With some extra tuning of n and e it is possible to improve on the computa¬ 
tional cost of the PHMC algorithm and a gain over the HMC algorithm can 
be obtained that can reach about a factor of 2. In particular it seems that 
when going to larger physical volumes this gain tends to increase. Another 
-substantial- gain can be obtained from the PHMC algorithm on massively 
parallel machines when several replica are run in parallel. 

3) Even if one decides to conservatively choose the polynomial parameters n 
and e, such that the computational cost becomes comparable to the one 
of the HMC algorithm, we still see a conceptual advantage of the PHMC 
algorithm. It samples configuration space differently from the HMC al¬ 
gorithm, allowing in particular for exceptionally small eigenvalues of the 
lattice Dirac operator to occur. Fermionic observables that are propor¬ 
tional to the inverse of these eigenvalues get corrected by the correction 
factor which makes the PHMC algorithm exact, yielding a finite contribu¬ 
tion to the (reweighted) sample average. We demonstrated this feature in 
a number of tests in this paper and showed that our way of treating these 
exceptional eigenvalues in the simulation is working in practise. If gauge 
conhgurations, carrying exceptionally small eigenvalues, are physically im¬ 
portant for some observables, the HMC algorithm, given its difficulty to 
generate such conhgurations, would have a very long autocorrelation time 
for these quantities. In this scenario the performance gain of the PHMC 
algorithm would be very large. Of course, an investigation of this issue is 
very expensive and should be performed -in our opinion- within projects 
aiming at the same time at physical results. 
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